% fig4a.m
% 演示在 (a, b) 参数平面内利用 parfor 并行扫描，
% 并对不同动力学状态（周期1~8、多周期MP、混沌CH、定态等）分配固定颜色。

clear; clc; close all;

%% 1. 固定参数（参考论文中设置）
alpha = 2.1;
beta  = 2;
gamma = 2;
Im    = 0.15;
F     = 0.15;

%% 2. 参数 (a, b) 扫描范围
Na = 200; Nb = 200;
a_range = linspace(0.57, 0.87, Na);
b_range = linspace(0.35, 0.65, Nb);

%% 3. 积分时间及去瞬态
T_total = 1000;
T_trans = T_total - 200;
optsODE = odeset('RelTol',1e-9,'AbsTol',1e-9);

% 扩展系统的初始条件：主系统 [q; sigma; phi] + 扰动向量 [v1;v2;v3]
X0 = [0; 0; 0; 1; 0; 0];

%% 4. 预分配存储状态标签的矩阵
stateMat = zeros(Na, Nb);

%% 5. 并行扫描 (a, b)
parfor iA = 1:Na
    ai = a_range(iA);
    localStateRow = zeros(1, Nb);  % 存储 (a固定,b变化) 的状态
    for iB = 1:Nb
        bi = b_range(iB);

        % 数值积分 扩展系统
        [t_sol, X_sol] = ode45(@(t,X) cnn(t, X, ai, bi, ...
            alpha, beta, gamma, Im, F), [0, T_total], X0, optsODE);

        % 提取主系统分量并计算 Vx
        q     = X_sol(:,1);
        sigma = X_sol(:,2);
        Vx = (alpha + beta*sigma).*q;

        % 只保留稳态段 (t >= T_trans)
        idx = (t_sol >= T_trans);
        Vx_ss = Vx(idx);

        % 用 findpeaks 寻找局部峰值
        [pks, ~] = findpeaks(Vx_ss);
        
        if isempty(pks)
            % 若没检测到峰值，认为定态
            localStateRow(iB) = 0; 
        else
            % 去重（微小数值波动可能导致虚假峰）
            pks_uniq = unique(round(pks, 3));
            nPeaks = length(pks_uniq);
            
            % 当 nPeaks==1 -> P1; nPeaks==2 -> P2; ... 
            % 超过8个唯一峰则认为混沌
            if nPeaks == 1
                localStateRow(iB) = 1;  % P1
            elseif nPeaks == 2
                localStateRow(iB) = 2;  % P2
            elseif nPeaks == 3
                localStateRow(iB) = 3;  % P3
            elseif nPeaks == 4
                localStateRow(iB) = 4;  % P4
            elseif nPeaks == 5
                localStateRow(iB) = 5;  % P5
            elseif nPeaks == 6
                localStateRow(iB) = 6;  % P6
            elseif nPeaks == 7
                localStateRow(iB) = 7;  % P7
            elseif nPeaks == 8
                localStateRow(iB) = 8;  % P8
            elseif nPeaks == 9
                localStateRow(iB) = 9;  % 多周期 (MP)，自行定义
            else
                localStateRow(iB) = 10; % 混沌 (CH)
            end
        end
    end
    stateMat(iA,:) = localStateRow;
end

%% 6. 定义离散 colormap，并与标签对应
myMap = [
    0.0  0.0  0.0;   % 0 -> black, 定态
    1.0  0.0  0.0;   % 1 -> red,   P1
    0.0  1.0  0.0;   % 2 -> green, P2
    0.0  0.0  1.0;   % 3 -> blue,  P3
    0.5  0.5  0.5;   % 4 -> grey,  P4
    1.0  1.0  0.0;   % 5 -> yellow,P5
    1.0  1.0  1.0;   % 6 -> white, P6
    1.0  0.0  1.0;   % 7 -> magenta, P7
    0.0  1.0  1.0;   % 8 -> cyan,  P8
    0.5  0.3  0.1;   % 9 -> brown, MP
    1.0  0.0  0.5;   % 10-> pink/red, CH
];

%% 7. 作图 (imagesc) 并使用离散 colormap
figure;
imagesc(a_range, b_range, stateMat'); 
set(gca, 'YDir','normal');  % 保证 b 从下到上增加
xlabel('a'); ylabel('b');
title('(a,b) 状态图');
clim([0, 10]);            % 状态标签 0~10
colormap(myMap);

% 设置 colorbar 的刻度与标签（需与 0~10 对应）
hc = colorbar('Ticks', 0:10, ...
    'TickLabels', {'ST','P1','P2','P3','P4','P5','P6','P7','P8','MP','CH'});

% 自定义颜色块高度，可结合 colorbar properties 调整
% set(hc, 'Limits', [0,10]) 等

